// STL
#include <algorithm>                              // std::max_element()
#include <iostream>                               // std::cerr
#include <fstream>                                // std::ofstream, std::ios
#include <functional>                             // std::function<>
#include <stdexcept>                              // std::runtime_error()
#include <string>                                 // std::string
#include <tuple>                                  // std::tie()
#include <vector>                                 // std::vector<>

// Boost
#include <boost/filesystem.hpp>                   // boost::filesystem::create_directories()

// jScience
#include "jutility.hpp"                           // sort_indices_desc()

// stats++
#include "statsxx/statistics.hpp"                 // statistics::mean()

// this
#include "statsxx/optimization/EA/EAParam.hpp"    // EAParam
#include "statsxx/optimization/EA/evolution.hpp"  // reproduce(), get_elite_individuals()
#include "statsxx/optimization/EA/fitness.hpp"    // calc_fitness()
#include "statsxx/optimization/EA/Individual.hpp" // Individual
#include "statsxx/optimization/EA/init.hpp"       // init_population()
#include "statsxx/optimization/EA/Population.hpp" // Population


//
// DESC: Perform an evolutionary algorithm.
//
// =====
//
// NOTE: See:
//
//          http://en.wikipedia.org/wiki/Evolutionary_algorithm
//
// =====
//
// TODO: NOTE: convergence() (subroutine) could probably be removed here in EA, replaced with convergence criteria.
//
// TODO: NOTE: This subroutine is (currently) not (even) used at all.
//
// TODO: NOTE: (The algorithm simply goes for the maximum number of generations.)
//
// -----
//
// TODO: NOTE: Perhaps should update to return multiple solutions.
//
// =====
//
// TODO: Some amount of carrying around duplicate information (e.g., in a Population and f[]) is being done; this should be updated.
//
inline std::vector<double> EA(
                              const std::function<std::vector<double>()>                            &new_individual,
                              // -----
                              const std::function<double(const int, const std::vector<double> &)>   &new_gene,
                              // =====
                              const std::function<double(const int, const std::vector<double> &)>   &mutate_gene,
                              // =====
                              const std::function<double(const std::vector<double> &)>              &fitness,
                              // -----
                              const std::vector<std::function<double(const std::vector<double> &)>> &inequality_constraints,
                              const std::vector<std::function<double(const std::vector<double> &)>> &equality_constraints,
                              // =====
                              const std::function<bool(const std::vector<double> &)>                &convergence_function,
                              // *****
                              const EAParam                                                         &param
                              )
{
    //=========================================================
    // ERROR CHECK
    //=========================================================

    if( !equality_constraints.empty() )
    {
        throw std::runtime_error("error in EA(): !equality_constraints.empty() --- and handling of these are not implemented!");
    }

    //=========================================================
    // INITIALIZE
    //=========================================================

    boost::filesystem::create_directories( param.out_dir );

    std::ofstream fitness_file((param.out_dir + "fitness.dat"), std::ios::trunc);

    //=========================================================
    // STEP 1. GENERATE INITIAL POPULATION
    //=========================================================

    Population population = init_population(
                                            new_individual,
                                            param
                                            );

    //=========================================================
    // STEP 2. CALCULATE INITIAL FITNESS
    //=========================================================

    std::vector<double> f;
    std::vector<double> v;
    std::tie(
             f,
             v
             ) = calc_fitness(
                              population,
                              // =====
                              fitness,
                              // -----
                              inequality_constraints,
                              equality_constraints
                              );

    std::vector<double> f_out;

    for( int j = 0; j < v.size(); ++j )
    {
        if( v[j] == 0. )
        {
            f_out.push_back(f[j]);
        }
    }

    if( !f_out.empty() )
    {
        fitness_file << "0   " << f_out.size() << "   " << statistics::mean(f_out) << "   " << *std::max_element(f_out.begin(), f_out.end()) << '\n';
    }
    else
    {
        fitness_file << "0   0" << '\n';
    }

    //=========================================================
    // STEP 3. EVOLVE
    //=========================================================

    for( int i = 0; i < param.max_gen; ++i )
    {
        //---------------------------------------------------------
        // STEP 3a & 3b. RE-PRODUCTION & RE-PRODUCE
        //---------------------------------------------------------

        Population offspring = reproduce(
                                         population,
                                         // -----
                                         f,
                                         v,
                                         // =====
                                         new_gene,
                                         mutate_gene,
                                         // =====
                                         param
                                         );

        //---------------------------------------------------------
        // STEP 3c. EVALUATE FITNESS OF NEW INDIVIDUALS
        //---------------------------------------------------------

        std::tie(
                 f,
                 v
                 ) = calc_fitness(
                                  offspring,
                                  // =====
                                  fitness,
                                  // -----
                                  inequality_constraints,
                                  equality_constraints
                                  );

        f_out.clear();

        for( int j = 0; j < v.size(); ++j )
        {
            if( v[j] == 0. )
            {
                f_out.push_back(f[j]);
            }
        }

        if( !f_out.empty() )
        {
            fitness_file << (i+1) << "   " << f_out.size() << "   " << statistics::mean(f_out) << "   " << *std::max_element(f_out.begin(), f_out.end()) << '\n';
        }
        else
        {
            fitness_file << (i+1) << "   0" << '\n';
        }

        //---------------------------------------------------------
        // STEP 3d. REPLACE INDIVIDUALS
        //---------------------------------------------------------

        population = offspring;

        // NOTE: The fitness values f[] (and constraint violations v[]) calculated (above) in Step 3c now correspond to "population".

        //---------------------------------------------------------
        // STEP 3e. CHECK FOR CONVERGENCE
        //---------------------------------------------------------

//        std::vector<Individual> champ;
//        get_elite_individuals(1, population, champ);
//
//        if( convergence_function(champ[0].genome) )
//        {
//            break;
//        }
    }

    //=========================================================
    // GET THE MOST FIT (FEASIBLE) INDIVIDUAL
    //=========================================================

    int n = 1;
    // -----
    bool feasible = true;

    std::vector<Individual> elite = get_elite_individuals(
                                                          n,
                                                          // -----
                                                          feasible,
                                                          // =====
                                                          population,
                                                          // -----
                                                          f,
                                                          v
                                                          );

    //=========================================================
    // CLEANUP & RETURN
    //=========================================================

    fitness_file.close();

    // RETURN
    if( !elite.empty() )
    {
        return elite[0].genome;
    }
    else
    {
        // TODO: NOTE: Perhaps a throw would be better here, so that the calling program could use try{}...catch{}?
        //
        std::cerr << "no feasible individual found in EA(): returning blank std::vector<>" << std::endl;

        return std::vector<double>();
    }
}
